In [1]:
import warnings
warnings.simplefilter("ignore")
#
import conda
conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
import os
#
import time
import datetime
import statsmodels.api as sm
import statsmodels.tools.eval_measures as smte
from cmocean import cm as cmo
import xarray as xr
import paramiko
import numpy as np
import pandas as pd
from netCDF4 import num2date
from netCDF4 import Dataset as NetCDFFile
#
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mc
import matplotlib.dates as mdates
from matplotlib.ticker import AutoMinorLocator
import matplotlib.path as mpath
from mpl_toolkits.basemap import Basemap
%matplotlib inline
#
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point
import cartopy.feature as feature
#
from cdo import *
cdo=Cdo()
In [2]:
plt.rcParams['figure.dpi'] = 300
plt.rcParams.update({'font.family': 'sans-serif', 'text.usetex': False,'pcolor.shading':'auto'})
In [3]:
if os.path.exists("format.py") == False:
os.system("ln -s ../../../anl/sample/format.py")
from format import *
mons=Mon_Name().Mon()
Before Computing EOFs¶
In [4]:
tmax=30*12
lonmin=0
lonmax=360
latmin=-90
latmax=-20
latmax2=-60
nlons=128
nlats=25
nlats2=18
numeof=10
infile="z500"
Making monthly climatology¶
In [5]:
if os.path.exists(infile) == False:
os.system("ngted -e AITM3:c'SFC1' ../cat/1dy/z500")
os.system("ngtavr ../cat/1dy/z500 -m 1mo -o ../cat/z500")
os.system("ngted -e AITM3:c'SFC1' ../cat/z500")
os.system("ln -s ../cat/z500")
if os.path.exists(infile+"monmean.nc") == False:
os.system("ngtcf "+infile)
cdo.settunits ('-days -settaxis,1981-01-15,12:00,1mon',
input=infile+".nc",output="tmp.nc")
os.system("mv tmp.nc "+infile+".nc")
cdo.ymonmean(input=infile+".nc",
output=infile+"monmean.nc")
cdo.settunits ('-days -settaxis,1981-01-15,12:00,1mon',
input=infile+"monmean.nc",output="tmp.nc")
os.system("mv tmp.nc "+infile+"monmean.nc")
Making anomaly from climatology¶
In [6]:
if os.path.exists(infile+"_S90-S20_anom.nc") == False:
cdo.sellonlatbox(str(lonmin)+","+str(lonmax)+","\
+str(latmin)+","+str(latmax),\
input=infile+".nc",\
output=infile+"_S90-S20.nc")
cdo.sellonlatbox(str(lonmin)+","+str(lonmax)+","\
+str(latmin)+","+str(latmax),\
input=infile+"monmean.nc",\
output=infile+"monmean_S90-S20.nc")
cdo.sub(input=infile+"_S90-S20.nc "+infile+"monmean_S90-S20.nc", output=infile+"_S90-S20_anom.nc")
Check trend¶
In [7]:
if os.path.exists(infile+"_S90-S20_anom_fldmean.nc") == False:
cdo.fldmean(input=infile+"_S90-S20.nc",output=infile+"_S90-S20_fldmean.nc")
cdo.fldmean(input=infile+"_S90-S20_anom.nc",output=infile+"_S90-S20_anom_fldmean.nc")
#
fldmean = NetCDFFile(infile+"_S90-S20_fldmean.nc" ).variables["zg"][:,0,0]
fldmeananom = NetCDFFile(infile+"_S90-S20_anom_fldmean.nc").variables["zg"][:,0,0]
time = np.arange(0,tmax,1)
#
plt.rcParams["figure.dpi"] = 300
fig = plt.figure(figsize=(6,2))
ax1 = fig.add_subplot(2,1,1)
ax2 = ax1.twinx()
ax1.plot( time,fldmean,
color="black",linewidth=1, label="raw")
ax2.plot( time,fldmeananom,
color="red",linewidth=1, label="anom")
ax2.axhline(y=0,color="k",ls="-",linewidth=0.5)
plt.xlim([1,tmax])
plt.xlabel("Number of mons")
plt.ylabel("z500")
#plt.legend(loc="best")
plt.show()
Check fields¶
In [8]:
def LatLon(nlons):
dlats = NetCDFFile(infile+"_S90-S20.nc").variables["lat"]
dlons = NetCDFFile(infile+"_S90-S20.nc").variables["lon"]
dlons2=np.zeros(nlons+1)
dlons2[0:nlons]=dlons[:]
dlons2[nlons]=-1*dlons[0]
flon,flat=np.meshgrid(dlons2,dlats)
return flon, flat
def LonLat2(nlons,nlats):
delta = 2.*np.pi/(nlons-1)
lons = (delta*np.indices((nlats,nlons))[1,:,:])
lats = -1*(0.5*np.pi - delta*np.indices((nlats,nlons))[0,:,:])
return lons,lats
def xy(lons,lats,m):
x, y = m(lons*180./np.pi, lats*180./np.pi)
return x,y
In [9]:
if os.path.exists(infile+"timeselmean.nc") == False:
cdo.timselmean(tmax,
input=infile+"_S90-S20.nc",
output=infile+"_S90-S20_timeselmean.nc")
mean=NetCDFFile(infile+"_S90-S20_timeselmean.nc").variables["zg"][0,:,:]
mean2=np.zeros((nlats,nlons+1))
mean2[:,0:nlons]=mean[:,:]
mean2[:,nlons]=mean[:,0]
if os.path.exists(infile+"timeselmean.nc") == False:
cdo.timselstd(tmax,
input=infile+"_S90-S20_anom.nc",
output=infile+"_S90-S20_timeselstd.nc")
std=NetCDFFile(infile+"_S90-S20_timeselstd.nc").variables["zg"][0,:,:]
std2=np.zeros((nlats,nlons+1))
std2[:,0:nlons]=std[:,:]
std2[:,nlons]=std[:,0]
monmean=NetCDFFile(infile+"monmean_S90-S20.nc").variables["zg"]
monanom=np.zeros((12,25,128))
anommean=np.zeros((25,128))
for mon in range(12):
monanom[mon,:,:]=monmean[mon,:,:]-mean[:,:]
monanom2=np.zeros((12,nlats,nlons+1))
monanom2[:,:,0:nlons]=monanom[:,:,:]
monanom2[:,:,nlons]=monanom[:,:,0]
Mean & standard deviation¶
In [10]:
plt.rcParams['figure.dpi'] = 300
fig = plt.figure(figsize=(6,4))
flon,flat=LatLon(nlons)
ax = plt.subplot(1,2,1)
ax.set_title("mean", fontsize=16)
m = Basemap(projection='spstere', lon_0=180,boundinglat=latmax)
x,y=m(flon,flat)
m.drawcoastlines()
cs = m.contourf(x,y,mean2)
cbar=plt.colorbar(cs, orientation='horizontal')
cbar.ax.tick_params(labelsize=6)
ax = plt.subplot(1,2,2)
ax.set_title("std", fontsize=16)
m = Basemap(projection='spstere', lon_0=180,boundinglat=latmax)
m.drawcoastlines()
cs = m.contourf(x,y,std2)
cbar=plt.colorbar(cs, orientation='horizontal')
cbar.ax.tick_params(labelsize=6)
fig.show()
Monthly anomaly¶
In [11]:
plt.rcParams['figure.dpi'] = 300
fig = plt.figure(figsize=(6,6))
lons,lats = LonLat2(nlons,nlats)
levels=np.arange(-200,220,20)
for i in range(12):
ax = plt.subplot(3,4,i+1)
ax.set_title(mons["name"][i], fontsize=16)
m = Basemap(projection='spstere', lon_0=180,boundinglat=latmax)
m.drawcoastlines()
cs = m.contourf(x,y,monanom2[i,:,:],cmap="coolwarm",levels=levels)
cax = fig.add_axes([0.1, -0.01, 0.8, 0.03])
art = plt.colorbar(cs, cax, orientation='horizontal')
art.ax.tick_params(labelsize=16)
plt.tight_layout(pad=0.2, w_pad=0.2, h_pad=0.2)
plt.show()
Compute EOFs¶
In [12]:
if os.path.exists(infile+"eofmonorg_pc_00000.nc") == False:
# compute modes
cdo.eof(str(numeof), \
input="-setmisstoc,0 "+infile+"_S90-S20_anom.nc", \
output=infile+"evalmonorg.nc "+infile+"eofmonorg.nc")
# compute principal coefficients
cdo.eofcoeff(input=infile+"eofmonorg.nc "+infile+"_S90-S20_anom.nc"\
, output=infile+"eofmonorg_pc_")
In [13]:
if os.path.exists(infile+"conmonorg.nc") == False:
#--read output data
ff = cdo.seltimestep("1/"+str(numeof), input=infile+"eofmonorg.nc", \
options='-f nc', returnCdf=True)
# read eigen vectors
eof = ff.variables["zg"][:,:,:]
# for variance calculation
cdo.remapcon("r"+str(nlons)+"x"+str(nlats),\
input=infile+"_S90-S20_anom.nc",\
output=infile+"conmonorg.nc")
#
if os.path.exists(infile+"evalmonorg_noweight.nc") == False:
if os.path.exists(infile+"variancemonorg.txt") == True:
os.system("rm "+infile+"variancemonorg.txt")
#
# Making a script
os.system("echo export CDO_WEIGHT_MODE=off > "
+infile+"variancemonorg.sh")
os.system("echo numeof=10 >> "
+infile+"variancemonorg.sh")
os.system("echo eval_file="+infile+"evalmonorg_noweight.nc >> "
+infile+"variancemonorg.sh")
os.system("echo eof_file="+infile+"eofmonorg_noweight.nc >> "
+infile+"variancemonorg.sh")
os.system("echo inp_file="+infile+"conmonorg.nc >> "
+infile+"variancemonorg.sh")
os.system("echo 'cdo eof,${numeof} ${inp_file} ${eval_file} ${eof_file}' >> "
+infile+"variancemonorg.sh")
os.system("echo 'cdo fldsum -timvar ${inp_file} tot1.nc' >> "
+infile+"variancemonorg.sh")
os.system("echo 'for ((i=1; i<=${numeof}; ++i)) ; do' >> "
+infile+"variancemonorg.sh")
os.system("echo 'cdo -seltimestep,${i} ${eval_file} tot2.nc' >> "
+infile+"variancemonorg.sh")
os.system("echo 'mode=`cdo output -mulc,100 -div tot2.nc tot1.nc | sed \"s- --g\"`' >> "
+infile+"variancemonorg.sh")
os.system("echo 'echo EOF MODE ${i}: ${mode} >>"+infile+"variancemonorg.txt' >> "
+infile+"variancemonorg.sh")
os.system("echo done >> "
+infile+"variancemonorg.sh")
# Run script
os.system("chmod 755 "+infile+"variancemonorg.sh")
os.system("./"+infile+"variancemonorg.sh")
os.system("rm tot1.nc tot2.nc")
In [14]:
cor=np.zeros(numeof+1)
jj=0
for line in open(infile+'variancemonorg.txt', 'r'):
jj=jj+1
va, vb, vc, vd= line[:-1].split(' ')
# print(va,vb,vc,vd)
cor[jj]=round(float(vd),2)
In [15]:
numeof2=3
labels=["SAM","PSA1","PSA2"]
In [39]:
sign=[1,1,-1,-1,1,1,-1,-1,1,1] # sign for each mode #
ff = cdo.seltimestep("1/"+str(numeof), input=infile+"eofmonorg.nc",\
options='-f nc', returnCdf=True)
# read eigen vectors
eof = ff.variables["zg"][:,:,:] # read eigen vectors
eof2=np.zeros((numeof,nlats,nlons+1))
eof2[:,:,0:nlons]=eof[:,:,:]
eof2[:,:,nlons]=eof[:,:,0]
In [40]:
pc2=np.zeros((tmax,numeof2))
detrend2=np.zeros((tmax,numeof2))
trend2=np.zeros((tmax,numeof2))
for ii in range(numeof2):
if os.path.exists(infile+"eofmonorg_pc_%5.5i_trend.nc" % (ii)) == False:
cdo.detrend(input=infile+"eofmonorg_pc_%5.5i.nc" % (ii),output=infile+"eofmonorg_pc_%5.5i_detrend.nc" % (ii))
pc = NetCDFFile(infile+"eofmonorg_pc_%5.5i.nc" % (ii)).variables["zg"][:,0,0]
pc2[:,ii] = (pc-pc.mean())/pc.std()*sign[ii]
detrend = NetCDFFile(infile+"eofmonorg_pc_%5.5i_detrend.nc" % (ii)).variables["zg"][:,0,0]
detrend2[:,ii] = detrend/detrend.std()*sign[ii]
trend2[:,ii] = pc2[:,ii] - detrend2[:,ii]
eof2[ii,:,:] = eof2[ii,:,:]*pc.std()
EOFs¶
In [31]:
def Color_Set(levels,cmapname):
import matplotlib.pyplot as plt
import matplotlib.colors as mc
cmap = plt.get_cmap(cmapname, len(levels)+1)
colors = cmap(range(cmap.N))
cmap, norm = mc.from_levels_and_colors(levels,colors,extend="both")
return cmap, norm
In [32]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import ListedColormap, BoundaryNorm, rgb2hex
def make_colorpalet(colorname="viridis",levels=[10,30,60,70,90]):
# https://hackmd.io/@h2tg95D2RP2ed-D8u-49Mg/S1moBqaRr
white="#ffffff"
lightcyan="#E0FFFF"
paleturquoise="#AFEEEE"
# color_code=[paleturquoise,lightcyan,white]
color_code=[paleturquoise,white]
num = len(levels)-len(color_code)+1
color_code.extend(get_color_code(colorname,num)[:])
# カラーマップ作成(→ 1)
cmap, norm = mc.from_levels_and_colors(levels,color_code,extend="both")
return cmap, norm
def get_color_code(cname,num):
cmap = cm.get_cmap(cname,num)
code_list =[]
for i in range(cmap.N):
rgb = cmap(i)[:3]
# print(rgb2hex(rgb))
code_list.append(rgb2hex(rgb))
return code_list
In [33]:
z500_levels = [-100, -70,-50,-40,-30, -20, -10, -5,-1,1,5, 10,20,30,40,50,70,100]
z500_cmapname = "twilight_shifted"
z500_cmap, z500_norm = Color_Set(z500_levels, z500_cmapname)
color_code=(get_color_code(z500_cmapname,len(z500_levels)+1)[:])
mid=(len(color_code)-1)/2
color_code[int(mid)]="#ffffff"
z500_cmap, z500_norm = mc.from_levels_and_colors(z500_levels,color_code,extend="both")
z500_label = "Geopotential height (500 hPa) [m]"
In [41]:
fig = plt.figure(figsize=(24,7))
lons,lats = LonLat2(nlons,nlats)
jj=0
for ii in range(numeof2):
#----draw eigen vector
jj=jj+1
ax = plt.subplot(1,numeof2,jj)
m = Basemap(projection='spstere', lon_0=180,
boundinglat=latmax)
m.drawcoastlines()
cs = m.contourf(x,y,eof2[ii,:,:],
norm=z500_norm,cmap=z500_cmap,levels=z500_levels)
ax.text(0.015,0.9, labels[ii] +
" ("+str(cor[ii+1])+"%)" ,
transform=ax.transAxes, fontsize=24,
bbox=dict(boxstyle='square',fc='w',
alpha=1.0), zorder=100)
cax = fig.add_axes([0.1, -0.05, 0.8, 0.05])
art = plt.colorbar(cs, cax,orientation='horizontal',ticks=z500_levels)
art.set_label(z500_label, fontsize=32)
art.ax.tick_params(labelsize=24)
plt.tight_layout(pad=0.2, w_pad=0.2, h_pad=0.3)
plt.show()
In [42]:
fig = plt.figure(figsize=(24,30))
jj=0
for ii in range(numeof2):
#----plot PC
jj=jj+1
ax = plt.subplot(numeof2,1,jj) #
f3 = cdo.selname("zg",input=infile+"eofmonorg_pc_%5.5i.nc" % (ii),
options='-f nc', returnCdf=True)
xx = np.arange(0,len(pc2[:,ii]),1) #
lns3 = ax.plot(xx, pc2[:,ii], "-", ms=10, color="black", linewidth=2)
ax.bar(xx,pc2[:,ii],width=1)
ax.plot(xx,trend2[:,ii],color="red",linewidth=3)
#
ax.tick_params(labelsize=24)
ax.set_xticks(np.arange(0,360,12))
ax.set_xticklabels(np.arange(1981,2011),rotation=70)
ax.set_ylabel(labels[ii]+"-index",fontsize=32)
ax.set_ylim((-3.0,3.0))
ax.axvline(x=120,ls='--',linewidth=1, color='k')
ax.axvline(x=240,ls='--',linewidth=1, color='k')
ax.axhline(y=1,ls='--',linewidth=1, color='k')
ax.axhline(y=0,ls='--',linewidth=1, color='k')
ax.axhline(y=-1,ls='--',linewidth=1, color='k')
plt.tight_layout(pad=0.2, w_pad=0.2, h_pad=0.5)
plt.show()
In [45]:
# 日付情報
time = NetCDFFile(infile+".nc").variables["time"][:]
date_offset = datetime.datetime(1981,1,1)
date_offset.day
y = np.zeros(tmax)
m = np.zeros(tmax)
d = np.zeros(tmax)
for tt in range(tmax):
time_data = date_offset+datetime.timedelta(days=time[tt])
y[tt] = time_data.year
m[tt] = time_data.month
d[tt] = time_data.day
np.savez("sammon.npz",sam=pc2[:,0], y=y, m=m, d=d)
np.savez("psa1mon.npz",sam=pc2[:,1], y=y, m=m, d=d)
np.savez("psa2mon.npz",sam=pc2[:,2], y=y, m=m, d=d)